^ im r. Ik;:! . mM\ 


NASA Technical Memorandum 102847 


€ 


A One-Equation Turbulence Transport 
Model for High Reynolds Number 
Wall-Bounded Flows 

Barrett S. Baldwin and Timothy J. Barth 


August 1 990 


(NASA-TM-1G2&47) A 

TRANSPORT MOU p L F nR 
wall-bounded FLOWS 


jN £- £ OU A T ION TUOBULE NC E 
i-iIGH REYNOLDS NUMBER 
(NASA) 23 p CSCL 200 


03/34 


NR1-10252 


Unci 3S 
03 10q31 


NASA 

National Aeronautics and 
Space Administration 


NASA Technical Memorandum 102847 


A One-Equation Turbulence Transport 
Model for High Reynolds Number 
Wall-Bounded Flows 

Barrett S. Baldwin and Timothy J. Barth, Ames Research Center, Moffett Field, California 


August 1 990 


NASA 

National Aeronautics and 
Space Administration 


Ames Research Center 

Moffett Field, California 94035-1000 



SUMMARY 

A one-equation turbulence model that avoids the need for an algebraic length scale 
is derived from a simplified form of the standard k — e model equations. After calibration 
based on well established properties of the flow over a flat plate, predictions of several 
other flows are compared with experiment. The preliminary results presented indicate that 
the model has predictive and numerical properties of sufficient interest to merit further 
investigation and refinement. The one-equation model is also analyzed numerically and 
robust solution methods are presented. 


INTRODUCTION 

One motivation for the developments documented in this report was the inability 
of well-established Navier-Stokes solvers using algebraic turbulence models to adequately 
predict several of the turbulent flow fields contained in the Viscous Transonic Airfoil Work- 
shopf. These flows contained significant separated flow and uniformly poor results were 
reported by all participants using the algebraic turbulence models of Baldwin-Lomax (ref. 
1) or Cebeci-Smith (ref. 2). (The results predicted by the ARC2D code were reported 
by Maksymiuk and Pulliam in ref. 3.) Two cases were adequately predicted by only one 
turbulence model, that of Johnson and King (ref. 4) as reported by King (ref. 5) and 
verified by Coakley (ref. 6) using an independent code. 

Another motivation is the need to treat flow problems in which multiple shear layers 
are present such that the determination of algebraic length scales is cumbersome and un- 
reliable. An example is the Coanda airfoil configuration computed by Pulliam, Jespersen, 
and Barth (ref. 7) which exploits tangential surface blowing. The need to avoid algebraic 
length scales leads to a consideration of k — e or related two-equation models. 

From our limited experience with two-equation models and reports by others (e.g., 
Sugavanum (ref. 8)) it became apparent that it would be worthwhile to investigate the 
possibility of transforming to variables other than the basic physical variables in an effort 
to avoid the well-known numerical difficulties that occur in the solution of the standard 
k — c equations. In the course of that investigation, a self-consistent one-equation model 
was found that also avoids the need for algebraic length scales. The main purpose of this 
report is to present the one-equation model and show the applicability of the model to a 
range of difficult turbulent flows. Results from computations of the two troublesome cases 
in the Viscous Transonic Airfoil Workshop are reported and significant improvement is 
achieved. 

The following two sections explain the rationale behind the development of the one- 
equation model. The next contains the solution for a self-similar turbulent wake to demon- 
strate a degree of generality of the model. In the overall design of the model equation, em- 
phasis was placed on numerical considerations so that extremely robust numerical solution 
methods could be used. The Numerical Implementation section gives some of the numer- 
ical theory and analysis needed to properly incorporate the one-equation model into flow 
solvers. The one-equation model has been implemented in a number of central-difference 

f Held in conjunction with the AIAA 25th Aerospace Sciences Meeting (January 1987). 
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and upwind finite-volume Navier-Stokes solvers in both two and three space dimensions 
and generalized coordinates. The computer code for these implementations can be ob- 
tained by contacting the second author (barth@prandtl.nas.nasa.gov). Finally in the 
appendix, we completely summarize the one-equation model for compressible flow. 

The authors are grateful to Drs. P.R. Spalart and T.J. Coakley for useful discussions 
and for reviewing the report. 

DERIVATION OF THE k - R T MODEL 


We begin with a standard form of the k — e equations (see Patel et al. (ref. 9)): 


Dk 

~Dt 


= V-(i/+— )Vk + P-e 


Dt 

m 




vt_ 


) V 6 + C tl 




( 1 ) 


where denotes the substantive derivative, + V ■ V and P the production 

of fc, P — u t It* - — f v t ■ From these two equations we are free to 

form a third by considering linear and nonlinear combinations. In our case, we do this to 


form equivalent systems which have improved numerical properties. We will return to this 
discussion in a later section. In particular, we consider a field equation for the “turbulence 
Reynolds number,” Rt 


Rt = — (turbulence Reynolds number) (2) 

ve 

The Rt field equation is obtained from the A: — e equations by considering differentials 
of J2t> dRr/Rr = 2 dk/k — de/e . It should be clear that the substantive derivatives 
of k and e as well as their respective source terms transform without approximation. In 
transforming the diffusion terms (which are modeled in both the k and e equations), we 
omit certain terms arising from the transformation to obtain a new diffusion model for the 
Rt equation. 

= (2 - c„)^P + (c„ - 2 )k + (v + — )V 2 (W2t) - i(V*) ■ V(vRt) (3) 

Lft le <T t <T € 

where 

v t = c^vRt) 

Note that since Ut should not depend on v at large Rt , the appropriate field variable is 
uRt = k 2 /e rather than Rt- 

Equation (2) can be rearranged in the form 

k a (k 1 + k 2 r 

uRt uRt 
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Without loss of generality, we can assign 


k\ — uRtP (at large Rt) (5) 

In regions where << fci this will result in P « e. Note that the relation 

e = < k , + h? = t?(i + = v r t p( i + ^) 2 

fcj All 

or 

k = y/uRrP (1 + ^) (6) 

is still completely general. Substitution of equation (6) into equation (3) and rearrangement 
leads to 

D ( U * T ) =(c e , - c tl WvR T P + (v + -)V 2 (vR t ) - — (Vi/t) • V(^r) 

i^r <r e cr € ^ 

~( 2 ~ c«t)r x i - y/vRrP- (2 - c ea )fc 2 

«1 + «2 

The systei 1 can be closed by substituting equation (6) in the k equation to obtain 


Dk 

~Dt 



k\ 

vRt 


+ V-(v + 


— )Vfc 


(8) 


k 2 = k — yf\ vR'j'P . 

Use of equations (4) and (5) requires further comment. These steps can be viewed as a 
means for exploiting the approximation resulting from equating production to dissipation, 
but in a way that rational procedures for departures from that approximation remain 
available. 

THE ONE-EQUATION R T MODEL FOR WALL-BOUNDED FLOWS 

It is of interest to note that by neglecting the last two terms in equation (7), a 
self-consistent one-equation model is obtained that should be a valid approximation over a 
major portion of shear layers. The main objective of this report is to develop modifications 
to this one-equation model that will allow it to be used in all parts of a shear layer. For 
that purpose, ail of the previous relations are taken to be applicable at sufficiently large 
Rt • To arrive at a model that is applicable in near-wail regions, the turbulence Reynolds 
number Rt is split into two factors 


r t = RTh(Rr) (9) 

where is a damping function such that Rt » Rt at large Rt . In addition, damping 
functions commonly used in k — e models (ref. 9) are introduced so that 


v t = vc^^Rt = VC^fpfsRT 


( 10 ) 
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and 


(ii) 


e' = e - D = 


(fci + fcj) 2 


vRt vf 3 R T 

The definition of ki applicable at small Rt (and at all Rt) is taken to be 

kj = vRtP 


( 12 ) 


This will allow k\ to be the dominant part of k in the near-wall region if the damping 
function fz is designed to accomplish that purpose. The resulting field equation for Rt is 


D(vRt) 


^ v 5_ D I i l/t \T7 t l{ ^ 


m = (Cej2 - C c x )y vRtP + (u + —)V 2 (uRt) - — (Vi/t) • V{uRt) (13) 
Ut (T € (t € 

Incompressible flow over a flat plate with zero pressure gradient is used to deter- 
mine suitable damping functions and to help calibrate the model. The thin shear layer 
approximation is used so that the production P is approximated by 


P = v t (u y ) 2 (thin shear layer assumption) 


(14) 


Then 


uR t P = c^uRruy) 2 fpf 3 
and the model equation reduces to 

DRt 


7-w ( c €a/2 c €i)y/ c ^f\ifzRT u y + (v + ){^T)yy ( z/ t)y(^T)y (15) 

Ut (T e <T € 

At sufficiently high momentum thickness Reynolds number the beginning of the log 
layer occurs in a region where the total shear stress is approximately constant and equal 
to the shear stress at the wall. In the log region and below, where advective terms are 
negligible, the x-momentum equation becomes simply 


1 


( v + v t )u y = u 2 (log region and below) 


(16) 


where u T is the friction velocity y/r wa u( p wa \u I** the log region where u v = u r /(«y), 
v << i/ t , and Vi = vc^Rt we have 


K 

v t = — = Ku r y 
U V 


where k is the Karman constant and 

Vt 


Rt = 


vc L 


k u T y 

Cp V 


(log region) 


(log region) 


In this case (i?T)yy = 0 and substitution in equation (15) with damping functions set to 
unity and advection terms zero produces the following formula: 


- = (* 
<?€ 




( 17 ) 
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In the region below the log layer an additional relation is needed to determine the 
definitions of Rt and } 3 {Rt)‘ After a study of the consequences, we have imposed 

Rt = — — — (log region and below) (18) 

Cft v 

We can ensure that this relation is consistent with the field equation for Rt (eqn. (15)) by 
requiring that the damping function / 2 (y + ) be adjusted to accomplish that purpose. By a 
procedure that will be described, the following damping functions have been determined: 

u =d 3 /d 3 

fz =D 2 D 3 

where 

£>i = l — exp(— y + /A + ), A + = 26 

D 2 = 1 - exp(-y + /A%), A\ - 10 

D 3 = 1 + B 3 exp(-y + /Aj)[ 1 - exp(-y + f A %) ], B 3 = 5.2, Aj = 15 

where y + = u T y/u . From these expressions we have that / M / 3 = D\D 3 . It should be 
emphasized that neither nor f 3 is equal to D\ or Z? 2 . However, it will be seen that 
implementation of the one-equation model requires only the evaluation of the product f tl f 3 
and does not require knowledge of / M and f 3 individually. Since D\ and D 2 are simpler 
functions than and f 3 , it is convenient to replace the product / M / 3 with D\D 2 wherever 
it occurs. 

Substitution of equations (17) and (18) in equation (15) with DRt/ D t = 0 yields 


/q 

(c ej / 2 - c (x )y/c^D^Dlyu y = (c t2 - c €l )^-f(u t )y 

K> 

Substitution of equation (16) to remove u y and equations (10) and (18) to remove v t 
produces 


M«+) = — +(1 - — K-t: + D l D 2 )(s/D7n 

Ce 2 C fJ Ky ' ' 


y 




(j+ ex P(~y + / A+ ) D 2 + -j+exp(-y + /A+) Di)) 


(19) 


For small y + we have the following limit value: 


lim / 2 (y + ) = ^ + (1 - 

y+— ►o C €3 c Ca 


« y/A+At 


After assignment of the parameters c tl ,c £j ,/c, equation (19) can be used to determine the 
damping function / 2 (y + ) when the combination 7^/3 = D\D 2 is known as a function of 
y"*". For our choice of constants we have a limit value of / 2 (0) ~ 0.781 . 
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From equation (18) it is seen that y + in the damping functions can be replaced by 


y 


+ = c -±r t 

K 


(when p x = 0) 


( 20 ) 


so that the damping functions can be expressed in terms of Rt or (implicitly) in terms of 
the turbulence Reynolds number Rt = Rt/ 3 {Rt )* For general flows that involve nonzero 
pressure gradients, the dependence on Rt or is not interchangeable. Experience will 
determine whether the functional dependence on Rt or is preferable. 

In this report, for simplicity, we have adhered to dependence on but we do not 
recommend this as a final choice. We have calibrated these functions by comparison with 
results from the (inner) Cebeci-Smith model (ref. 2) using equation (16). The resulting 
prediction of u t from that model in the log region and below is 

(vt)cs = ^ [\/l + (4 «J/ + Di) 2 - l] 

For the present model, substitution of equation (20) in equation (10) yields 

Vt = v k y + = V DiD 2 K y + 

For either model u + is obtained by the integration of equation (16): 


u 


+ 



dy + 

1 + u t fv 


Figure 1 shows a comparison of the predicted variations of v t jv and u + . The close agree- 
ment of the two models is not surprising because the damping functions D\ and D 2 were 
designed for that purpose. 

From equations (10), (12), (14), (16), and (18) explicit formulas for kf and P + in the 
inner region are obtained 


1 y/PiZ?2 k y + 
y/c^ 1 + I>1 D 2 Ky+ 


(log region and below) 


P + = 


Ky^PiPi 

(1 + K y+ Dj D 2 y 


(log region and below) 


A plot of the kf variation is shown in Figure 2. It is somewhat surprising that ki provides a 
more realistic turbulent energy distribution than several k—e models, although no attention 
was paid to that purpose in the design of D\ and D 2 . 


6 



Figure 1, Comparison of U + and 
Vt between Cebeci- Smith and the one- 
equation model. 



Figure 2. Turbulent Kinetic Energy 
and damping function (f 2 ) distribu- 
tions. 


Since has been introduced, it seems necessary that it be fully defined. The above D$ 
function vas calibrated to produce a reasonable variation of e in the near-wall region, 
although neither e nor D$ is needed in calculations based on our one-equation model. To 
complete the definition of Rt and the quantity D in equation (11) is taken to be uk yy 
(from Mansour, Kim, and Moin (ref. 10)). 

For the sake of completeness, plots of RtRti /3> and e are shown in figures 3-5, 
although, as mentioned earlier, they are not needed for implementation of our model. The 
e distribution was evaluated assuming that k « ki, (i.e., k 2 << ki) in the log region and 
below. The damping function D$ was designed to make the e distribution resemble that 
found from the direct simulations reported in reference 10. However the resulting rather 
realistic distribution is a byproduct that was not anticipated. 




y 


+ 


Figure 3. Rt , Rt distributions. Figure 4. Damping functions. 
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Figure 5. Near- wall epsilon distribution. 


To complete the calibration of the one-equation model, values of the parameters 
k, c^, c tl , c €2 are needed. In this report we have adopted the widely used values k = 
0.41, Cp = 0.09 and set c €j = 2.0, so that the last term in equation (7) will be identically 
zero. We have adjusted c €l to match the calculated skin friction coefficient for incompress- 
ible flow over a flat plate to a compilation of experimental data that is well represented by 
the Karman-Schoenherr formula given in Hopkins and Inouye (ref, 11): 


~r = 17.08(log 10 Reef + 25.11 log 10 Reg + 6.012 (Karman — Schoenherr formula) 

C f 


The value of c £l currently used is 1.2, which is well below the “standard” value of 1.44 but 
is not the lowest on record, see reference 9. An alternative to our low value of c £l would be 
to impose a relatively large value of Rt in the outer flow. From equation 13 it can be seen 
that a constant value of Rt is self-preserving in regions where production is zero (unlike 
k and e in the k — e equations). However, in this report, we adhere to c Cl = 1.2 and a low 
level of Rt in the free stream such that Ut in the free stream (and at the leading edge of 
an airfoil) will not be large compared to the molecular viscosity. 
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Figure 0 . Law-of-the-wall. Figure 7. Skin friction on flat plate. 

In this report, we extend the model to compressible flow following a practice which 
has been used for algebraic models (ref. 2). We assume that all the above relations apply 
to compressible flow with a temperature-dependent u and with 


Pt = pvt 


Figure 6 graphs computed law-of-the-wall solutions (tt + = u* ju T versus y + ) for subsonic 
and supersonic Mach numbers where u* is defined (see, for example, Rubesin and Horstman 
(ref. 12)): 

fU 

u* = I y/pjp^ai dU 

Jo 

For incompressible flow this plot corresponds to the conventional law of the wall. The paper 
by Hopkins and Inouye (ref. 11) contains several alternative procedures for applying the 
Karman-Schoenherr formula to compressible flow. In the procedure by Sommer and Short 
an empirically determined temperature T aa is defined: 

T aa = T edge ( 1 + 0.035M e 2 dae ) + 0A5{T waH - T edge ) 

The momentum thickness Reynolds number Re$ is then adjusted to a corresponding in- 
compressible value Re$: 

Re 0 = Re 0 p(T edge )/p(T aa ) 

Finally, the resulting C / from the Karman-Schoenherr formula is adjusted to the com- 
pressible value by multiplying by the temperature ratio T edge /T aa . Figure 7 compares skin 
friction for compressible flow over a flat plate with the foregoing “theory” at subsonic and 
supersonic Mach numbers. The agreement for both subsonic and supersonic Mach num- 
bers is good. The excellent agreement at low Mach numbers was expected because the 
model was calibrated in the incompressible limit. 

Figure 8 contains law-of-the-wall plots for flat plate boundary-layer flow computed 
using several values of free-stream R? • Note that vtjv = c^Rt = 0. 09 Rt when damping 
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functions are_at unity. The plotsjndicate insensitivity to free-stream values for all values of 
free stream Rt tested, except (Rt)o o = 100 which is the only value at which v t jv exceeds 
unity in the free stream. 



Figure 8. 


Sensitivity of boundary layer to values of free-stream Rt- 


To further assess the performance of the one-equation model, computations were per- 
formed for two of the troublesome cases in the Viscous Transonic Airfoil Workshop. Figures 
9-11 show the mesh and solution for viscous flow over the RAE 2822 airfoil geometry at 
Moo = 0.75, a — 2.72°, and a Reynolds number of 6.2 million. Computed lift and drag 
coefficients for the one-equation model were Cl — 0.771, Co — 0.0352. This is compared 
with the Baldwin-Lomax solution (also plotted) which produced lift and drag coefficients 
of Cl = 0.895, Co = 0.0279. The shock position is substantially improved and the overall 
pressure coefficient distribution is in much better agreement with experiment. 



Figure 9. RAE 2822 mesh. 



0.75). 
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Figure 11. Pressure coefficient comparison. 


The i.ext geometry and flow conditions provide a much more severe test case for 
turbulence modeling. The geometry is the standard NACA 0012 computed at — 0.799, 
a = 2.26°, and a Reynolds number of 9 million. Figures 12-14 show the grid and solution 
at these flow conditions. Figure 14 plots C p distributions for the one-equation model 
(Cl — 0.340, Co = 0.035), the Baldwin-Lomax model (Cl = 0.531, Co — 0.048), and the 
one-equation model with advection removed from the model equation on upper-surface 
flow only (Cl = 0.589, Cp = 0.048). 



Figure 13. Mach contours (M a 0 = 
Figure 12. NACA 0012 mesh. 0.8). 
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Figure 14. Pressure coefficient comparison. 

The improvement in the full one-equation model is dramatic. The upper-surface shock 
wave is moved forward almost 20 percent chord and is in good agreement with experi- 
ment. The solution obtained with and without upper-surface advection clearly indicates 
the importance of upstream influences obtained via advective terms in this separated flow. 
Note that we have used the corrected angle of attack suggested by the experimenters 
(a exp == 2.86°). The discrepancy in lower-surface pressure and the effect of upper-surface 
shock location seems to indicate that this angle-of-attack correction does not adequately 
account for wind tunnel wall interference. 

THE ONE-EQUATION R T MODEL FOR FREE SHEAR LAYERS 


In this section we examine the one-equation Rt model for the self-similar turbulent 
wake flow. Following the order-of-magnitude arguments of Tennekes and Lumley (ref. 13), 
the x momentum equation for an eddy viscosity model calculation reduces to the following 
simplified form: 

UU X = (u t Uy)y (21) 

Using similar arguments, equation (7) reduces to 


t)x =(c € 2 — c €l ) 


^ + i m r t )„ - m,(Rt),] 

v a € 





where u t = vc^Rt, P = ^t(U y ) 2 , and fcj = yj 

We are interested in the one-equation model that results from neglecting the last 
two terms of equation (22). However, it is worthwhile to gain an appreciation of how 
that approximation can apply more generally than in regions where production equals 
dissipation. Wake flows provide an interesting example. In this flow, production is zero 
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at the centerline, but dissipation is not. Moving away from the centerline, the assumption 
that production is equal to dissipation becomes valid. Fortunately, even at the centerline 
of a wake flow we can rationalize neglecting the last two terms in equation (22). Clearly at 
the centerline of a wake the first of these two terms is identically zero because production 
vanishes there. Also note that many k — e modelers choose c c 2 = 2, which would remove 
the last term identically. The primary question still remains as to the overall validity of 
the one-equation approximation in a simple wake. Assuming self similarity of a turbulent 
wake, we can analyze the situation in detail. Following the scaling procedure of Tennekes 
and Lumley for wake flows, we set 

U = Uo + U 8 f(t) 


r t = ^Mf) 


(23) 

(24) 

vt = vc^Rt = U,lcnh(£) (25) 

where £ = y/l, U s = Ax _1 / 2 , l = Bx 1/2 , and U s « U 0 are assumed. Substituting in 
equations (21) and (22) (with the last two terms of eqn. (22) omitted) we obtain 

W + */«) + <*(*/«)« = 0 (26) 

Hht + ^ [hh U - (h<) 2 ] + -\U\h = 0 (27) 

( T € 

where 


0 


U 0 B 

2A 


.08 


(28) 


The numerical value of [3 was evaluated from experiment by Tennekes and Lumley, based 
on the definitions U s = \U — f7o|moi» an d by choosing a length l as the distance from the 
centerline such that | U — Uo\/\U — Uo \ max — exp(-l/2). 

If h ( o is a known function, equation (26) can be evaluated in the form 


/ = exp 



(29) 


This expression is then differentiated to replace \f{\h with (3£f /c^. Equation (27) can then 
be expressed in the following form: 


‘.-m 


(30) 


Pi(0 


Qi(0 


where fi j = 02=0 (t 5- ) • With Pi and Qi considered known, the resulting equation 

is linear in and the solution is 


-s: 


Qi{Oexp[Pu(t) - Pii(t)} d Z 
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where Pu = Pi By iterating the quadratures for / and h, a numerical solution 
can be obtained in about 10 steps. In figure 15 we plot the solution of this equation and 
compare it with results from the constant eddy viscosity solution as well as results from 
experiment. 
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Figure 15. Comparison of experiment, and theory for self-similar wake flow. 

The data points are from Townsend (ref. 14). The dashed curve was computed using 
our one-equation model with the values of c € j and c € 2 from the original “standard” k — e 
model. The solid curve is based on the values of parameters that provide the best results 
for the wall-bounded flows considered in previous sections. The departure of the solid 
curve from the data is no worse than that of the dotted curve, which is from the widely 
used (in wakes) constant eddy viscosity approximation. 

NUMERICAL IMPLEMENTATION 

In this section we consider discretization of the one-equation model. We begin in a 
rather abstract setting by introducing the notion of positive operators for discrete systems. 
This will allow us to construct extremely robust algorithms for the one-equation model. 
Current implementations of the one-equation model in two- and three-space dimensions 
employ an implicit factored ADI solver for the scalar equation (decoupled from the mean 
flow equations). In the following analysis, implicit unfactored schemes will be considered. 
Good numerical behavior of this system is a prerequisite for the more general situation 
which might include matrix factorization. 

We first define the solution vector TZ on a two-dimensional, logically rectangular mesh 
where 72^ & jfcj) with 

TZ = [7^1,1,72.1,2, ...,7 £1 ,m, 7£ 2| i > 7£ 2> 2, •••,7£ 2i m, ■•• ? 7£;v,i,72jv,2, —,TZn,m ] T • 

As we will see, our one-equation model with discretized advection and diffusion terms 
produces a system of ordinary differential equations of the form 



where M('R, ) is a matrix operator (possibly nonlinear) representing the discretization of 
advection and diffusion and D is a diagonal matrix with positive entries representing the 
source term. We can construct implicit and explicit schemes of the form 

.71 . Tl-\- 1 yTl ^Tl yTl 

[i + At0M(n )}{k -n ) = a t(-M(n ) + d)-r ( 32 ) 

or after rearrangement 

.n . n-(- 1 y 

[I+At9M(TZ )}n =[I -{\-6)AtM{n ) + AtD]K ( 33 ) 


for all 9 € [0,1]. As we will show, we can design numerical approximations for advection 



and diffusion which guarantee the following properties of ): 


1 . M('R-) is a diagonally dominant monotone (M-type) matrix. M-type matrices are diag- 
onally dominant matrices with positive diagonal entries and negative off-diagonal entries. 


2. M(H) has zero row sum. This will be due to the use of the chain rule form of the 
equations 


A well-known property of M-type matrices is that they have non-negative inverses, i.e., 
elements of the inverse are non-negative. We see from property 1 that the left-hand-side 
matrix of (33) 



is an M-type matrix and consequently 


[I + A i6M(n n )}- 1 >0, 9 > 0 

The right-hand-side matrix of (33) 


[I - (1 - 0)AtM{n n ) + AtD] 


(34) 


is also unconditionally non-negative for 6 = 1 and non-negative under a CFL-like condition 

,o 

for 1 > 6 > 0. Assuming 'R, >0, from the solution update equation 

.71 4-1 vTl y 71 y n 

n =[I+At6M(n )]-*[/ - (1 - 6)AtM{K ) + A tD]U 
,1 1+1 

positivity of TZ is guaranteed whenever (34) is a nonnegative operator. 

We now investigate stability properties of the numerical scheme. Given properties 1 
and 2, it is a simple matter (see, for example, Barth and Lomax (ref. 15)) to show that 


||[/+Af0M(^V 1 || oo <l 
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We also have that 


||[J-(l-0)A*M(ft)]||oo<l 

under the CFL-like condition for positivity. If we first consider the system without the 
source term ( D — 0) and Dirichlet boundary conditions, we have stability under the CFL- 
like condition for non-negativity 

^ 71 + l y 71 . , n 

\\n ||oo =\\[i + At 8 M(n )}- 1 [i-(i-$)AtM(Ti)}n h*. 

<||[/ + A^M(^ n )]- 1 |U ||[/-(l-<?)A/M(^)]|| 00 ||^ n || 00 

<ll^ n |loo 

In the presence of the linear source term with D > 0, we obtain the following stability 
estimate: 


^n +1 ^ n n 

II ||oo<(l+At||I>|| 0O )||ft II oo < (1 + Atmajc(D J:J ))|| K ||oo 

3 

This result is expected because the differential equation admits growth of this sort in the 
presence of the source term (with positive coefficient). 

We now turn to the actual discretization of the individual terms. For convenience we 
rewrite the diffusion terms (assuming v is constant) 

(u + — )V 2 (uRt) - — (Vi/t) • V(vRt) = 2(u + — )V 2 (vR t ) - — V • v t V(t ;R T ) 

<T t <T ( <T e (T ( 


Using this identity, we rewrite equation (13) ( S 2 = 

+ V.VH t ) = 2( 1 / + -)V 2 Ht)--V^ ( V{vR t ) + (Ce, h - c £l )^JJ~ 3 SR t 

OX (T € (T € 

We first approximate the advective terms by using a standard, first-order accurate upwind 
approximation (u ± = (u ± |u|)/2,u ± = (t? ± \v\)/2): 


V • V(uR t ) ttaZllj+i'k + 0ZHj,k + llKj- 1,* 


where 


a? = u - t , 

a Ax 


7a = ft = ~( a « + 7a) 




A />'•*’ 


7 a = 


vt„, « = -K+7i!) 


Ay ’■ k ’ 

Note that this discretization automatically satisfies both properties 1 and 2 mentioned 
previously. 
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The diffusive/antidiffusive terms are approximated by central differencing: 


2 (v + V 2 (uR t ) » 2 (*/ + . 


Tlj+ijk — + R^j- i,fc 

Ax 2 

'fcj.fc+i ~ + 'Rj.fe-i 

Ay 2 


V ■ — V(i/J 2 t ) 

<7 




Ax 2 


+ 


(gjj+l ~ ^j,fc) ~ 1/2 ^j.fc-i) 

Ay 2 


By combining these two expressions we obtain 

2 [v + -) v 2 (^t) - V • — v(uRt) ~<* d n j+ltk + pin j%k + 7^-1,* 

\ <7 / <T 

+ a d^j,fc+i + + Td^j,fc-i 


where 


Ax 2 

2 (^ + -) -(-) 

V cr / \<T/j+i/ 2 ,k 

75 = Ax 2 

[2 (*/+-) -(-) 

L v <r/j,fc \<T / j-\/2,k 

y 1 

Ay 2 

2(l /+-) -(-) 

V cr/j,fc v < r/ i,fc+l/2. 

ld ~ Ay 2 

2 { V+ ~) “(“) 

V cr t j,k . v a /j,k-i/2\ 


« — («s + t 5 ) 

« = -(«S+-r 5 ) 


Note that in this form we do not have automatic satisfaction of property 1 which would 
require that both a and 7 be positive coefficients (property 2 is automatically satisfied). 
It is important to realize that if these coefficients become negative this effect is entirely a 
result of poor grid resolution. To see this, assume a smooth variation of u t and expand v t 
in a Taylor series. For example we have that 


1/2, k = + O(Ax) 


and 


= 


Ay 2 


2 u+(~) + 0 { Ax) 

\<r / j,k 


From this equation it is clear that this coefficient is guaranteed positive for small enough 
Ax. If we assume that ut >> u, we obtain the following restriction for non-negativity of 
a and 7: 

(^t)j+l,fc < ^( U t)j,kl (ut)j-l t k 3 (vt)j,k 

(vt)j,k+ 1 < 3(i/ t )j,fc, (vt)j,k- 1 < 3 (^t)i.fc 
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In our implementation of the algorithm we strictly enforce this condition. Whenever these 
conditions are violated we limit the amount of anti-diffusion added so as to maintain 
positivity. Keep in mind that the better remedy, however, is usually grid refinement. 

We conclude this section with a final remark concerning grid-resolution requirements 
for the one-equation model. Because the variable Rt was designed to behave linearly in 
the near-wall region for zero pressure gradient boundary layers (as previously shown in 
Figure 3 and equation (18)), we typically only require a mesh wall spacing comparable 
to the Baldwin-Lomax model (y+ < 3.5). This limit is required for accurately estimating 
wall shear which is used in the present implementation to determine in the damping 
functions. This is a drastic improvement over near- wall formulations of the k — e equations 
which typically require mesh wall spacing less than y + = 0.2 as mentioned in reference 9. 
This removes from the flow solver much of the stiffness resulting from the extremely fine 
meshes which must be used in the near- wall formulations of the k — e model. In figure 16 we 
plot the law of the wall for the flat plate using several mesh wall spacings 0.5 < y+ a u < 3.1 
to demonstrate the solution independence with wall spacing. 



Figure 10. Sensitivity of boundary layer to wall spacing. 


CONCLUSIONS 

A one-equation turbulence model has been introduced that avoids the necessity for 
finding an algebraic length scale. The preliminary results presented indicate that the model 
has predictive and numerical properties of sufficient interest to merit further investigation 
and refinement. 
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APPENDIX - SUMMARY OF THE ONE-EQUATION MODEL 


This appendix gives a complete summary of the one-equation model. We begin with 
the field equation for Rt : 


E&k) = _ c ^n- + (l/ + ^)V ! („fi T ) - — (V^O • V(vRt ) (35) 

In this equation, we use the following functions: 

= ( C €2 ~ C €l)\/C~ii/ K 

<?€ 

vt =c^(uRt)D 1 D 2 
Pt =pvt 

Di =1 - exp (— y + /A + ) 

D 2 =1 - exp (-y + /A%) 


(9Ui 

our 

\ SOi 2 I 

OU k \ 2 

\0*i 

dx iy 

1 dxj 3 ‘ ' 

\dx k ) 


h(y + ) =— + (1 - — )(^T + D,D 2 )(y/D^Dl 

C (3 C f2 Ky+ 

eM ~ y+ ,A+) D% + ir xp(_!,+M2+) Di)) 

For all calculations we have used the following constants: 


k =0.41, c fl = 1.2, c €j = 2.0 
Cfl =0.09, A + = 26, A+ = 10 

We also recommend the following boundary conditions for (35): 

1. Solid Walls: Specify Rt = 0. 

2. Inflow (V ■ n < 0): Specify Rt — {Rt ) oo < 1. 

3. Outflow (V • n > 0): Extrapolate Rt from interior values. 
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